//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : misc_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef MISC_KERNELS_H
#define	MISC_KERNELS_H

namespace common
{

// kernel to initialize device array
template <class T>
__global__ void kernel_memset(T* dst, const int n, const T val)
{
    int tid = bx * Bx + tx;
    if (tid < n)
        dst[tid] = val;
}

// kernel to go from embedding grid to the r1/r2/b1/b2 storage
template <class T>
__global__ void kernel_from1to4(T *dst, const T *src, 
    const Grid grid1, int nx0)
{   
    int cgx = grid1.cx;
    int cgy = grid1.cy;
    int lda_new = grid1.nx;
    int lda_old = nx0;

    int u = BORDER_WIDTH + bx * DIM_COMPUTE_BLOCK + tx;
    int v = BORDER_WIDTH + by * DIM_COMPUTE_BLOCK + ty;
    int loc;

    for (int k = 0; k < T1TO4_BF; ++k)
    {
        T val = src[v * lda_old + u];

        loc = lda_new * BORDER_WIDTH
                    + lda_new * by * (DIM_COMPUTE_BLOCK >> 1)
                    + lda_new * (ty & 1) * (BORDER_WIDTH + cgy)
                    + lda_new * ((DIM_COMPUTE_BLOCK / T1TO4_BF) >> 1) * k
                    + lda_new * (ty >> 1)
                    + BORDER_WIDTH
                    + bx * (DIM_COMPUTE_BLOCK >> 1)
                    + (tx & 1) * (BORDER_WIDTH + cgx)
                    + (tx >> 1);

        dst[loc] = val;

        v += (DIM_COMPUTE_BLOCK / T1TO4_BF); // (v += By, By = blockDim.y)
    }
}

// kernel to go from the r1/r2/b1/b2 storage to the embedding grid
template <class T>
__global__ void kernel_from4to1(T* dst, const T* src, 
    const Grid grid1, const int nx0)
{
    int cgx = grid1.cx;
    int cgy = grid1.cy;
    int lda_new = grid1.nx;
    int lda_old = nx0;

    int u = BORDER_WIDTH + bx * DIM_COMPUTE_BLOCK + tx;
    int v = BORDER_WIDTH + by * DIM_COMPUTE_BLOCK + ty;
    int loc1;
    int loc2;

    for (int k = 0; k < T1TO4_BF; ++k) 
    {
        loc2 = lda_new * BORDER_WIDTH
                    + lda_new * by * (DIM_COMPUTE_BLOCK >> 1)
                    + lda_new * (ty & 1) * (BORDER_WIDTH + cgy)
                    + lda_new * ((DIM_COMPUTE_BLOCK / T1TO4_BF) >> 1) * k
                    + lda_new * (ty >> 1)
                    + BORDER_WIDTH
                    + bx * (DIM_COMPUTE_BLOCK >> 1)
                    + (tx & 1) * (BORDER_WIDTH + cgx)
                    + (tx >> 1);


        T val = src[loc2];

        loc1 = lda_old * v + u;

        dst[loc1] = val;

        v += (DIM_COMPUTE_BLOCK / T1TO4_BF);
    }
}

// kernel to divide r2 nodes from a grid over r1/r2/b1/b2 nodes of the
// next higher (coarser) grid (1 level down)
template <class T>
__global__ void kernel_divideR2in4(Grid grid1, const Grid grid0)
{   
    int lda_old = grid0.nx;
    int cgx_old = grid0.cx;
    int cgy_old = grid0.cy;
    int cgx     = grid1.cx;
    int cgy     = grid1.cy;
    int lda_new = grid1.nx;

    int u = BORDER_WIDTH2 + cgx_old + bx * DIM_COMPUTE_BLOCK + tx;
    int v = BORDER_WIDTH2 + cgy_old + by * DIM_COMPUTE_BLOCK + ty;
    int loc0;
    int loc1;

    for (int k = 0; k < R2TO4_BF; ++k) 
    {
        loc0 = lda_old * v + u;

        T val_CC = grid0.cc[loc0];
        T val_NN = grid0.nn[loc0];
        T val_EE = grid0.ee[loc0];
        T val_SS = grid0.ss[loc0];
        T val_WW = grid0.ww[loc0];
        T val_NE = grid0.ne[loc0];
        T val_SE = grid0.se[loc0];
        T val_SW = grid0.sw[loc0];
        T val_NW = grid0.nw[loc0];

        loc1 = lda_new * BORDER_WIDTH
                    + lda_new * by * (DIM_COMPUTE_BLOCK >> 1)
                    + lda_new * (ty & 1) * (BORDER_WIDTH + cgy)
                    + lda_new * ((DIM_COMPUTE_BLOCK / R2TO4_BF) >> 1) * k
                    + lda_new * (ty >> 1)
                    + BORDER_WIDTH
                    + bx * (DIM_COMPUTE_BLOCK >> 1)
                    + (tx & 1) * (BORDER_WIDTH + cgx)
                    + (tx >> 1);

        grid1.cc[loc1] = val_CC;
        grid1.nn[loc1] = val_NN;
        grid1.ee[loc1] = val_EE;
        grid1.ss[loc1] = val_SS;
        grid1.ww[loc1] = val_WW;
        grid1.ne[loc1] = val_NE;
        grid1.se[loc1] = val_SE;
        grid1.sw[loc1] = val_SW;
        grid1.nw[loc1] = val_NW;

        v += (DIM_COMPUTE_BLOCK / R2TO4_BF);
    }
}

// kernel to substitute the r1/r2/b1/b2 nodes from a grid into the r2 nodes 
// from the next lower (finer) grid (1 level up)
template <class T>
__global__ void kernel_compress4toR2(Grid grid1, const Grid grid0)
{
    int cgx_old = grid1.cx;
    int cgy_old = grid1.cy;
    int cgx     = grid0.cx;
    int cgy     = grid0.cy;
    int lda_old = grid1.nx;
    int lda_new = grid0.nx;

    int u = BORDER_WIDTH2 + cgx_old + bx * DIM_COMPUTE_BLOCK + tx;
    int v = BORDER_WIDTH2 + cgy_old + by * DIM_COMPUTE_BLOCK + ty;
    int loc1;
    int loc0;

    for (int k = 0; k < R2TO4_BF; ++k) 
    {
        loc0 = lda_new * BORDER_WIDTH
                    + lda_new * by * (DIM_COMPUTE_BLOCK >> 1)
                    + lda_new * (ty & 1) * (BORDER_WIDTH + cgy)
                    + lda_new * ((DIM_COMPUTE_BLOCK / R2TO4_BF) >> 1) * k
                    + lda_new * (ty >> 1)
                    + BORDER_WIDTH
                    + bx * (DIM_COMPUTE_BLOCK >> 1)
                    + (tx & 1) * (BORDER_WIDTH + cgx)
                    + (tx >> 1);

        T val_CC = grid0.cc[loc0];
        T val_NN = grid0.nn[loc0];
        T val_EE = grid0.ee[loc0];
        T val_SS = grid0.ss[loc0];
        T val_WW = grid0.ww[loc0];
        T val_NE = grid0.ne[loc0];
        T val_SE = grid0.se[loc0];
        T val_SW = grid0.sw[loc0];
        T val_NW = grid0.nw[loc0];

        loc1 = lda_old * v + u;

        grid1.cc[loc1] = val_CC;
        grid1.nn[loc1] = val_NN;
        grid1.ee[loc1] = val_EE;
        grid1.ss[loc1] = val_SS;
        grid1.ww[loc1] = val_WW;
        grid1.ne[loc1] = val_NE;
        grid1.se[loc1] = val_SE;
        grid1.sw[loc1] = val_SW;
        grid1.nw[loc1] = val_NW;

        v += (DIM_COMPUTE_BLOCK / R2TO4_BF);
    }
}

} // namespace common

#endif // MISC_KERNELS_H

